### CLUSTER ANALYSIS

# working directory -------------------------------------------------------
setwd("C://Users/hpeter.INTRANET/Desktop/multivariate stats R/2024/scripts/")


# load data ---------------------------------------------------------------
load("Doubs.RData")  


# Load required packages
library(vegan)


# Remove empty site 8
spe <- spe[-8, ]
env <- env[-8, ]
spa <- spa[-8, ]
latlong <- latlong[-8, ]


## Hierarchical agglomerative clustering of the species abundance 
## data

# Compute matrix of chord distances among sites
spe.norm <- decostand(spe, "normalize")
spe.ch <- vegdist(spe.norm, "euc")

image(as.matrix(spe.ch))
heatmap(as.matrix(spe.ch))
?heatmap


# Compute single linkage agglomerative clustering
?hclust
spe.ch.single <- hclust(spe.ch, method = "single")


# Plot the dendrogram using default options
par(mfrow=c(3,2), mar=c(4,4,1,1))
plot(spe.ch.single, main = "Chord - Single linkage", xlab="")


# Compute and plot complete-linkage agglomerative clustering
spe.ch.complete <- hclust(spe.ch, method = "complete")
plot(spe.ch.complete, main = "Chord - Complete linkage", xlab="")


# Compute UPGMA clustering
spe.ch.UPGMA <- hclust(spe.ch, method = "average")
plot(spe.ch.UPGMA, main = "Chord - UPGMA", xlab="")


# Compute UPGMC clustering
spe.ch.centroid <- hclust(spe.ch, method = "centroid")
plot(spe.ch.centroid, main = "Chord - Centroid", xlab="")


# Compute Ward's minimum variance clustering
spe.ch.ward <- hclust(spe.ch, method = "ward.D2")
plot(spe.ch.ward,  main = "Chord - Ward", xlab="")




#k-means clustering (non-hierarchical)
spe.ch.k<-  kmeans(spe.ch, centers=3)
spe.ch.k


# Cophenetic correlations -------------------------------------------------
# Single linkage clustering
spe.ch.single.coph <- cophenetic(spe.ch.single)
cor(spe.ch, spe.ch.single.coph)


# Complete linkage clustering
spe.ch.comp.coph <- cophenetic(spe.ch.complete)
cor(spe.ch, spe.ch.comp.coph)


# Average clustering
spe.ch.UPGMA.coph <- cophenetic(spe.ch.UPGMA)
cor(spe.ch, spe.ch.UPGMA.coph)


# Ward clustering
spe.ch.ward.coph <- cophenetic(spe.ch.ward)
cor(spe.ch, spe.ch.ward.coph)


# Shepard diagram ---------------------------------------------------------
par(mfrow=c(1,1), mar=c(4,5,3,1))
plot(spe.ch, spe.ch.single.coph,
  xlab = "Chord distance",
  ylab = "Cophenetic distance",
  asp = 1, xlim = c(0, sqrt(2)),
  ylim = c(0, sqrt(2)),
  main = c("Single linkage", paste("Cophenetic correlation =", round(cor(spe.ch, spe.ch.single.coph), 3))))
abline(0, 1)
lines(lowess(spe.ch, spe.ch.single.coph), col = "red", lwd=3)



## TASK 1: Use the oribatide mite dataset, calculate a similarity measure and perform clustering analysis.
# Q1: based on visual inspection of the dendrogram - which clustering algorithm produces the "best" results?
# Q2: how many groups of sites - based on mite community similarity are there? 
# Q3: explore the impact of log-transformation of mite abundances on the cluster analysis? 
#     What does this tell us in terms of the importance of abundant mite species for the overall structure of the mite community?







# find optimal number of clusters -----------------------------------------

library(NbClust)
?NbClust

plot(spe.ch.single)
Nb.single<-NbClust(spe, diss=spe.ch, distance = NULL, min.nc=2, max.nc=16, 
             method = "single", index = "ch")
Nb.single
plot(2:16,Nb.single$All.index, xlab="number of clusters", ylab="Calinski and Harabasz index")
abline(v=10, col="red", lty=2, lwd=2)




# advanced plotting -------------------------------------------------------
library(dendextend)

UPGMA.dend <- as.dendrogram(spe.ch.UPGMA)   #convert from hclust object to dendrogram
plot(UPGMA.dend)                            #plot dendrogram



Nb.UPGMA<-NbClust(spe, diss=spe.ch, distance = NULL, min.nc=2, max.nc=16, 
                   method = "average", index="ch")
Nb.UPGMA      #there are 2 optimal groups.

colors_to_use <- Nb.UPGMA$Best.partition      #define colors and sort according to tips in dendrogram
colors_to_use<-colors_to_use[order.dendrogram(UPGMA.dend)]

labels_colors(UPGMA.dend) <- colors_to_use   #change color of tip labels depending on best partition


#change color of branches using 5 (k-means) clusters
labels_colors(UPGMA.dend)<-1
UPGMA.dend <- UPGMA.dend %>% color_branches(k = 5)
plot(UPGMA.dend)


# TASK 2: Use NbClust to determine the optimal number of clusters in the mite dataset (see previous task) - how does this agree with your visual assessment?